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The optical response and the ground state of graphene and graphene-like systems are determined selfcon- 
sistently. Deriving equations of motion for the basic variables, graphene Bloch equations are introduced and 
combined with a variational Ansatz for the ground state. Within the Hartree-Fock approximation, this approach 
reproduces the gap equation for the ground state. The results show that the Coulomb interaction drastically 
influences the optical response of graphene and introduces an extremely sensitive dependency on the dielec- 
tric environment via screening. Regarding the effective fine-structure constant as control parameter, a transition 
from a semimetal to an excitonic insulator is predicted as soon as the effective graphene fine-structure constant 
exceeds a value of roughly 0.5. Above this critical value, the computed optical spectra exhibit a pseudogap and 
several bright p-like excitonic resonances. 
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I. INTRODUCTION 

The research interest in graphene has increased dramati- 
cally since its first isolation by Geim and co-workers in 2004. 
In particular its highly unusual electronic and optical proper- 
ties has led to a cascade of both theoretical and experimental 
investigations, exploiting its fundamental underlying physics 
as well as potential applications in electronic and optoelec- 
tronic devices. 

The key for understanding the unique electronic properties 
of graphene is its exotic band structure that differs substan- 
tially from most other condensed matter systems. Based on a 
tight-binding (TB) model, Wallace predicted 1947' an elec- 
tronic single-particle spectrum exhibiting two distinct cross- 
ing points. In the vicinity of these so called Dirac points, the 
dispersion is a cone similar to the light cone in relativistic me- 
chanics, with the Fermi velocity vp replacing the speed of 
light. The occurrence of the cones results from the symmetry 
between the two equivalent sublattices that build the honey- 
comb lattice. The sublattice wave functions can be combined 
into a pseudo-spinor which then obeys the ultra-relativistic 
Dirac equation. Hence, from a QED point of view, electronic 
excitations close to the Dirac points of graphene can be con- 
sidered as charged, massless, chiral fermions. From a con- 
densed matter point of view, as the density of states vanishes 
at the Dirac points, graphene can be considered either as a 
semimetal or a vanishing gap semiconductor 

However, as the tight binding Hamiltonian neglects many 
body interactions completely, the role of the electron-electron 
Coulomb interaction is still not well understood and the sub- 
ject of ongoing research.^"^^ As a convenient measure of the 
relative importance of the Coulomb interaction, one can use 
the effective fine-structure constant ac ~ e? jehvp, where e 
is the effective background dielectric constant. For freestand- 
ing graphene in vacuum olq 2.41, indicating prominent 
Coulomb interaction effects. 

Generally, one can distinguish between Coulomb modifica- 
tions of the electronic ground state and Coulombic signatures 
in the excitation properties. From strongly correlated systems 



it is known that the Coulomb interaction can induce a transi- 
tion from a semimetal to a Mott insulator where, unlike in con- 
ventional semiconductors, the gap results from the electron- 
electron rather than the electron-ion interactions. In particu- 
lar, the electron exchange interaction has been identified as the 
dominant mechanism responsible for the opening of a gap.^** 

Methodically, most theoretical treatments of the ground 
state of a many body systems either use highly simplified 
model Hamiltonians or they rely on either perturbative or vari- 
ational approaches. For graphene, perturbative studies based 
on a renormalization-group analysis predict a logarithmic di- 
vergence of the Fermi velocity, stabilizing the semimetallic 
ground state. However, nonperturbative methods yield a 
semimetal to insulator transition at sufficiently high coupling 
strengths*'"'" where predicted critical values for ac range from 
0.5 to 1.5. In particular, the possibility of an excitonic conden- 
sate has been explored. '^''^ 

From semiconductor physics, it is known that the excita- 
tion spectrum in the vicinity of the fundamental band gap 
is dominated by Coulomb bound electron-hole pairs leading 
to the appearance of excitonic resonances in the absorption 
spectra. Mathematically, the excitons obey a hydrogen-like 
Schrodinger equation, which is known as Wannier equation. 
The exciton binding energy in typical semiconductors is in 
the range of 1 to 100 meV, i.e. much smaller than the hy- 
drogen ground binding energy. However, bound states exist 
even in the presence of strong static (background) screening 
which preserves the long-range 1/r Coulomb tail. In metals, 
the highly mobile caiiiers effectively screen this long ranged 
part of the Coulomb interaction with the consequence that ex- 
citonic effects are of minor importance. In graphene, the ex- 
istence of bound pair states and the importance of screening 
effects are subjects of current research. Since real massless 
particles in nature are neutral, the Dirac two-body problem 
has become a topic of interest only within the graphene re- 
search and, until recently, it has not even been clear if bound 
Dirac pairs exist. '^•^'" 

Experimentally, an important tool to study microscopic pro- 
cesses in many-body systems is optical spectroscopy, provid- 
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FIG. 1. The lattice and the reciprocal lattice. 



ing information both on the system ground state and the exci- 
tation properties. Theoretically, several methods to model the 
optical response of a quantum mechanical many-body system 
are well established, e.g. the density matrix approach, non- 
equilibrium Green functions, or the systematic cluster expan- 
sion approach.^' In general, these methods rely on the knowl- 
edge of the initial state, which is usually the ground state. 

In this paper, we extend the microscopic approaches and 
develop a framework that allows us to determine the ground 
state and the optical response of graphene and graphene-like 
systems on the same level of approximation. Our method com- 
bines the equations of motion with a variational approach. 
Within the Hartree-Fock approximation, we obtain gap equa- 
tions that determine a Bogoliubov ground state.""''** 

With this initial state, we then calculate the linear optical 
response. The resulting spectra show several bright excitonic 
resonances below the pseudo-gap. 

The work is organized as follows. In Section II, we present 
the Hamiltonian within the tight-binding approximation. Us- 
ing this basis, we derive the equations of motion for the 
dynamical quantities of interest and discuss the role of the 
Coulomb interaction. Similar equations have been derived 
previously for carbon nanotubes^^'"'^. In section III, we apply 
the variational principle to the total Hamiltonian and impose 
the equations of motion as constraint to obtain a stationary 
ground state. The resulting gap equations are presented and 
analyzed numerically. In Section IV, we derive the graphene 
Bloch equations and compute the linear optical response. 



II. THE HAMILTONIAN AND EQUATIONS OF MOTION 

The system Hamiltonian for a single graphene sheet inter- 
acting with a classical light field can be written as 



H = Ho + Hi + Hc 



Here, 

Ra.Rb ^ 

describes the motion of the electrons of mass mo in the peri- 
odic lattice potential, V{x — R) is the effective core potential 
of the carbon atom located at R, and {Ra} and {Rb} are the 
coordinates of the carbon atoms on each sublattice (See Fig. 
1). 
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(2) 



describes the electron-electron interaction via the Coulomb 
potential V{x — x') — /e\x — x'\ where e is the dielectric 
constant of the environmental medium. 



Hr = / d'^a; 



2moc 



i;^x)^p-A{x) + A{x)-p--A'}'iP{x) (3) 



describes the light-matter interaction within the minimal cou- 
pling substitution and A is the vector potential for the optical 
field. 

The interaction Hamiltonian couples the dynamics for the 
vector potential to the expectation value of the particle cur- 
rent that can be calculated from the Heisenberg equation of 
motion. 
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for the relevant operator O, yielding a set of coupled differ- 
ential equations. As is well known, the two-particle Coulomb 
interaction couples the dynamics for the iV-particle expecta- 
tion values to those of the N + 1-particle expectation values, 
which is known as the hierarchy problem. In order to achieve 
a closed set of equations, this hierarchy must be truncated. A 
systematic truncation scheme is provided by the cluster ex- 
pansion, which has been proven to work quit well under many 
different excitation conditions.-" 

In principle, the system dynamics can be described by solv- 
ing the resulting coupled set of differential equations starting 
from a predetermined initial state. In a typical experimental 
setup, the system is excited from the ground state and the re- 
sponse to an externally applied field is measured. Thus, to 
analyze such experiments theoretically, it is crucial to have 
an adequate description of the ground state. In semiconductor 
physics, band structure calculations usually provide a suitable 
basis to expand the field operators and a good approximation 
for the ground state. However, this may not be the case in sys- 
tems where strong carrier-carrier Coulomb interactions influ- 
ence the ground state properties. 

In the following, we therefore follow an approach where 
we treat Hq within the tight-binding approximation and use 
the resulting eigenfunctions as basis to represent the total 
Hamiltonian. Within this basis, we then derive the equations 
of motion on the singlet level, which is equivalent to the 
time-dependent Hartree-Fock approximation. To determine 
the ground state on the same level of approximation, we apply 
the variational approach for the system energy and impose the 
equations of motion as constraints to guarantee stationarity of 
the ground state. 
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FIG. 2. Schematic plot of the yr-valence and the tt* -conduction band 
structure of graphene according to the tight-binding model 



A. Tight-binding Hamiltonian 

Following the tight-binding approach, we expand the field 
operators in terms of the carbon wave functions, 



ip{x) 



1 ^ e^''-R^<P{x-RA)ak 

k,RA 



— V e"'- 

k.Rs 



ijA{x) +ijB{x). (4) 



Here, (bk) annihilates a particle in the state {fc} on the 
sublattice A (B), and 4){r) = r e"''/^'' cos z9/\/327r(F is the 
carbon 2pz orbital responsible for the optical and electronic 
electronic properties of graphene. The parameter d — aB/Z^s 
controls the effective spreading of the carbon wave functions. 

Inserting Eq. (4) into the Hamiltonian and taking into ac- 
count nearest neighbor hopping only, one obtains for the 
single-particle part 



Hq — ^ Ep 

k 



d\dk 



lf{k)aibk + jr{k)bldk 



(5) 

Here, Epis the expectation value of Ho, 7 is the matrix ele- 
ment describing hopping processes between neighboring car- 
bon atoms. Si = f(l,V3), ^3 = f(l,-\/3), and ^3 = 
(—a, 0) are the vectors connecting any carbon atom to its three 
next neighbors (see Fig. 1), a is the nearest neighbor distance 
and 



/(fc) = E' 



ik-Si 



is a function that depends on the lattice symmetry properties 
only. Diagonalization of Hq yields the TB single-particle band 
structure 



e: 



EF±l\fik)\ 



(6) 



shown if Fig. 2. 

Within the tight binding approximation, the shape of the 
bands depends on the lattice symmetry only. The valence and 
conduction bands are symmetric and touch each other at the 
two nonequivalent points — (27r/3a, ±27r/3\/3a) of 

the Brillouin zone, where the function /(fc) — 0. Since the 



ground state of Eq. (5) has a completely filled valence and an 
empty conduction band, the nodes occur exactly at the Fermi 
level Ep. In the vicinity of these nodes, a first order Taylor 
expansion gives 



k) 



3a 



Thus, the dispersion is a cone with linear coefficient 37a/2 = 
hvp, similar to the light cone in relativistic mechanics. Con- 
sequently the electrons behave like massless Dirac fermions 
with the Fermi velocity replacing the speed of light. Both ex- 
perimental and theoretical data give a value of approximately 
10® m/s for the Fermi velocity. 

In the following, we will treat Hq within the TB approxi- 
mation. Within the band structure picture, the Hamiltonian is 
represented in terms of the electron and hole operators 



efc = (flfe + f{k)bk 
hU = ^{-dk + fik)bk), 



(7) 
(8) 



with /(fe) = /(fc)/|/(fc)|, yielding 
Ho^Y.^\f{k)\[eiek + hl^h^k] 

k 



7l/(fc)|). (9) 



Here, the second term only contributes if the electron and hole 
symmetry is broken, e.g. by doping, and the last term repre- 
sents the energy of the filled valence band. Since this term is 
constant, it can be omitted. Taking the tight-binding ground 
state as a reference, the energy with respect to this reference 
is obtained by normal ordering the Hamiltonian within the 
electron-hole picture. 

Important Coulomb contributions arise from scattering pro- 
cesses where each electron remains on its specific sublattice, 
defining the generic matrix element 



V{q) 



F{qd) 



that depends on the momentum transfer only and is given by 
the 2D bare Coulomb potential modified by the background 
dielectric constant e and a form factor F{qd). The form factor 



F{qd) 



dV 



(j3^'giq-(p-p')g-'3|2-2'l 



with F{Q) ~ 1 results from the finite extension of the carbonic 
p^-orbitals perpendicular to the plane and decreases monoton- 
ically with q. Hence, the finite value of d = as/Zctf can be 
interpreted as the effective thickness of the graphene sheet. 
As an intrinsic length scale, it fixes the graphene energy unit 
Eq = hvp/d and is crucial for obtaining finite values for the 
exciton binding energy. ""^ The Coulomb matrix elements for 
the processes where at least one electron is scattered from one 
sublattice to the other are much smaller and vanish exactly at 
the Dirac points. These will be neglected. 
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Expanding the field operators in terms of the electron and 
hole wave functions produces 2"^ different contributions, some 
of which describe equivalent processes. The 10 nonequivalent 
normally ordered contributions to the Coulomb interaction 

qhh' 

qkk' 
qhk' 

+ 5Z '^kk'^^)^i + q^-kh-k'+qek' 
qkk' 

+ \T. ^kk'i^)A + q^i'-qhlk'hU 
qkk' 

+ ^ X! ^kk'i^y^-k-qh-k'+qCk'ek 
qkk' 

- Y.^kk'{QA + qhlkA-^q^k 
qkk' 

qkk' 

+ 51 ^kk'i<l)f^-kh-k'+q'h-k-qek' 

qkk' 

+ E Vki'iQ)el.,hl^hl^,h^^^„ (10) 

qkk' 

have three different matrix elements 

X (l±r(fc'-q)/>')) (11) 

Vki'iq)^\v{q) (r(fc + q)/»-l) 

X (l + .f(fc'-g)/>')) (12) 

that not only depend on the momentum transfer but also on 
the momenta of the involved scattering particles. Only within 
the linear approximation, translational invariance is recovered 
and the matrix elements with q = k' k, relevant for the 
Hartree-Fock approximation, are given by 

V^{k - k') = \v{\k - k'\) (1 ± cos{9 - 9')) (13) 

^^(fc - k') = ±^y(|fc - k'\) sin(6l - e'). (14) 

The scattering processes characterized by the different ma- 
trix elements are schematically presented in Fig. 3. The re- 
pulsive e-e and h-h and the attractive e-h scattering where 
both particles remain in their initial bands are described by 
(Fig. 3a). The corresponding processes, where both parti- 
cles change their bands, are described by V~ (Fig. 3b). Both 
y+ and are symmetric with respect to fc — fc' and equal 
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FIG. 3. Schematic representation of the Coulomb matrix elements 



for both Dirac points. Auger processes, where an electron- 
hole pair is created or annihilated under simultaneous energy 
and momentum transfer to other quasiparticles, are not ruled 
out by energy and momentum conservation and described by 
the matrix elements (Fig. 3c). These matrix elements are 
antisymmetric with respect to k — k' and have opposite sign 
for the two distinct Dirac points. As usual, processes with zero 
momentum transfer describe a divergent self-interaction that 
cancels with the electron-ion and ion-ion interaction in the 
jellium limit and are explicitly subtracted from the Coulomb 
Hamiltonian. 

When computing the optical response, the light-matter in- 
teraction is often treated within the dipole approximation. In 
crystalline solids, the dipole approximation can be obtained 
by a multipole expansion and subsequent coarse graining over 
an elementary lattice cell, taking only monopole and dipole 
contributions into account. In typical direct-gap semiconduc- 
tors, the distinct symmetry properties of the valence and con- 
duction bands allow us to associate monopole contributions to 
intraband- and dipole contributions to interband transitions, 
respectively. As inter- and intraband transitions involve very 
different energy scales in a large-gap semiconductor, an opti- 
cal field couples only to interband transitions while intraband 
transitions are in the teraherz range. 

In graphene, the situation is quite different. As both the 
TT-valence and tt* -conduction band are constructed from pz 
atomic orbitals, they have the same angular momentum quan- 
tum numbers and monopole and dipole contributions do not 
distinguish between different bands. Hence, dipole transitions 
necessarily involve a superposition of fc-states required to 
build a p-like collective state. Additionally, due to the van- 
ishing gap at the Dirac points, inter- and intraband transitions 
take place on the same energy scale. 

To derive the correct interaction Hamiltonian, we start from 
Eq. (3), apply the Coulomb gauge 'V ■ A — 0, and expand the 
field operators in terms of the tight binding wave functions. 
Making the assumption that the vector potential varies slowly 
on the length scale of the lattice constant a and sheet thickness 
d, one finds 

k 

- —Y (n{k)-Aaibk+h.c.) (15) 



5 



f4 




B 



FIG. 4. Schematic representation of the Ught-matter interaction 



where A denotes the field at the position of the graphene 
sheet, z = 0. The first term in the first line of Eq. (15) de- 
scribes the interaction of the vector potential with the directed 
current hkpk/m(^, where = —e{a\ak+h\hk) is the charge 
density. This term vanishes because of parity. The second term 
in the first line of Eq. (15) is proportional to the total elec- 
tron density rio of the vr-band, yielding a Drude-like response. 
The remaining term describes true quantum mechanical tran- 
sitions, whose strength is determined by the optical matrix el- 
ement 

/2 A'f 

(16) 



and 



M 



(17) 



On-site momentum matrix elements vanish because of parity 
and hence, in graphene, optical transitions involve inter sub- 
lattice hopping processes. Within the linear approximation, 
the dipole matrix elements in the vicinity of the Dirac points 
are given by 



7r±(fe) 



'2M 



-i-n 1 6 



u± 



with u± ~ ± ie.y)/ \/2, showing that the two degenerate 
K points couple to circular polarization components of the 
optical field. 

Transforming the light-matter interaction Hamiltonian into 
the electron-hole picture yields within the linear approxima- 
tion. 



H 



[p] 



moc 



± i- 



.e TT 



cos e {elek + h^_Ji-k - l) 

k 

^A±sin0(4/il^-/i_feefc) (18) 



moC 



where the ± sign refers to the distinct Dirac points and circular 
polarization components respectively. Eq. (18) clearly shows 
that optical excitations not only involve interband, but also 
intraband transitions. Furthermore, the angle dependence in 
the light-matter Hamiltonian assures that only p-like states 
couple to an external optical field. 



B. Equations of Motion 

In this section, we use the basic single-paiticle expectation 
values 



fk = {e\ek) = {h^_Ji-k) 
Pk (h-kek) 



(19) 



as dynamical variables and derive their equations of motion 
within the time-dependent Hartree-Fock approximation. Due 
to the band symmetry, the electron and hole populations are 
equal. Furthermore, as the two equations of motion for the 
two distinct K points are related by the parity transformation 
k ^ —k, — > A^, we restrict our discussion to a single 
Dirac point. 

Evaluating the commutators in the Heisenberg equations of 
motion, we obtain the closed set of differential equations 



ih$-Pk = 2SfcPfc - (1 - 2fk)nk - ih^Pk 



At 



dt 



fi^jk = -2im[Pink]-hj^fk 

Here, we introduced the notation 



(20) 



(21) 



coll 



Efe = hvpk - J2 - fc') -V-ik- fe')] fk' 

k' 

- i V l/^(fc - k') ImPk' + ^ A± cos 9 
k^ moc 

= ^l[fk]+Ai:k[Pk,PLA] (22) 

rik-^Y.^^^'^- fc')^fc' + v-ik- k')p*, 

k' 

.elTrl 



T i^A^ sin0 - 2 y VCk - k')fk' 
moc "rf 

k 

= n§ [Pk,p;:,A]^iAnk [fk]- 



(23) 



In Eqs. (20) and (21), the terms d/dt|coii refer to incoherent 
scattering contributions beyond the Hartree-Fock approxima- 
tion. To solve the coupled equations, we have to supply the 
appropriate initial and boundary conditions. 

The generalized Rabi energy flk consists of two contribu- 
tions, il^ contains the external field and the internal field of 
the polarization Pk and differs from the standard Rabi fre- 
quency known from semiconductors only by the anisotropy 
resulting from the chiral nature of the Dirac electrons. The ad- 
ditional part Ail contains contributions from the populations 
fk via the Auger processes. Similarly, the generalized renor- 
malized energy has not only terms due to the populations 
contained in S'', but also Auger contributions from the po- 
larizations and an energy renormalization due to the external 
field in AS^. For large gap semiconductors, these contribu- 
tions can be neglected because of energy and momentum con- 
servation. However, for systems with small or vanishing gap, 
as in graphene, the basic quantities / and P are inherently 
mixed, due to the Auger terms. 
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Close to the K point, the Auger contributions to and E 
can be written as 

Arife = ± ^ V{k - k') sm{9 - e')h' 

k' 

AY.k = ± ^ V(fe - fc') sm{e - e') ImPk' , 

k' 

respectively. For isotropic densities / and P, the angle inte- 
gration will give zero, and the Auger contributions vanish. The 
light-matter terms, however, are in themselves anisotropic. 
The light will therefore excite angular-dependent densities for 
which the Auger contributions are finite. 

In semiconductor physics, the transverse optical field di- 
rectly couples only to the nondiagonal expectation values P 
and P*, which has the meaning of an optical polarization and 
can be probed by optical spectroscopy. Unlike in graphene, the 
diagonal populations fk affect optical semiconductor spec- 
tra only indirectly via their coupling to the polarization in 
the semiconductor Bloch equations (SBE).^** As the linear 
coupling to the optical field results in a radiative decay of 
any macroscopic polarization on a picosecond time scale,^"* 
a (meta-)stable initial state may be prepared to contain inco- 
herent diagonal populations but no macroscopic polarizations. 
In particular, choosing the initial fk — models excitation 
from the ground state. Expanding the solution of the SBE into 
powers of the exciting fields gives the well known series that 
contains only the odd powers for the polarization and the even 
powers of the optical field for the populations. Hence, for 
weak optical fields, the equations for the polarizations and 
populations can be solved iteratively. The homogeneous part 
of the polarization equation produces the generalized Wan- 
nier equation with the Pauli blocking prefactor (1 — 2/). This 
homogeneous part is responsible for the formation of bound 
excitons at low densities, that appear as resonances of the ex- 
cited system, energetically below the band gap. Inserting in- 
creasing amounts of incoherent initial populations leads to a 
gradual bleaching of excitonic resonances in the absorption 
spectra and finally produces gain once inversion is reached. 

A similar behavior might be expected for the case of 
graphene. However, the interpretation of the single-particle 
energy renormalization as "band edge" renormalization is 
problematically in a gapless system. Also, bound exciton so- 
lutions of the Wannier equation would predict resonances at 
negative energies, which is clearly unphysical. This patholog- 
ical behavior would show up for any incoherent initial pop- 
ulation that allows for bound-state solutions of the general- 
ized Wannier equation. Consequently, such a population can- 
not correspond to a physically meaningful initial state of the 
system. 

In particular, if the linear Wannier equation 



2hvFk(f>xk - [V+{k - k' 



k' 

Ex 0, 



'^\k' 



V-{k-k') 



'\k'. 



Xk ■ 



(24) 



has bound state solutions, the tight-binding ground state with 
fk = Pk = cannot correspond to the ground state of 
the Coulomb-interacting system. The linear Wannier equation 



has been analyzed in Ref. 30. Here, it has been shown that 
bound-state solutions indeed exist for an effective graphene 
fine-structure constant exceeding the critical value ac ~ 0.5. 
Below this value, the Coulomb interaction is too weak to cre- 
ate bound states. In this weakly interacting regime, Eqs. (20) 
and (21) can be solved directly to obtain the optical response 
and the tight-binding Dirac sea can serve as ground state of 
the system. In the regime with strong Coulomb interactions, 
Eqs. (20) and (21) are not the true graphene Bloch equations, 
since they do not describe the optical response defined with 
respect to the correct many-body ground state. 



III. GROUND STATE 

A. Derivation of the Gap Equations 

A convenient method to determine the ground state is pro- 
vided by the variational principle 



6{H) 



A=0 



= 0, 



(25) 



where S{H)\a=q is the expectation value of the total energy 
in the absence of external fields. Due to the Coulomb inter- 
action, the expectation value of the Hamiltonian contains the 
same many-body correlations as the Heisenberg equations of 
motion and cannot be solved exactly. To be consistent with 
the equations of motion, we express the energy expectation 
value in terms of the same variables as the equation of mo- 
tion. Within the mean field approximation, these are /fc, Pj, 
and P^, yielding the variational equation 

= ^ {2tkSfk ~ niSPk - ClkSP^) . (26) 

k 

where barred quantities refer to ground state expectation val- 
ues. However, the dynamic variables are coupled and thus can- 
not be varied independently of each other 

In order to be radiatively stable, the macroscopic current 
associated with the ground state must vanish. This condition 
is fulfilled for all isotropic distributions fk and P^. Due to 
this isotropy, the related Auger contributions to and Cl^ 
vanish. Furthermore, the ground state clearly has to be a sta- 
tionary state. Therefore, the expectation value of any operator 
in this state has to be static, such that the time derivative of the 
basic variables has to vanish. Thus, we impose the stationary 
solution of the equations of motion as constraints to assure a 
steady state. 



= 2EfcPfe - (1 - 2fk)nk 
= -2Im 



(27) 
(28) 



Equation (27) defines the possible distributions for the non- 
diagonal populations P^ in the presence of a given isotropic 
diagonal population fk subject to the condition that P^ is real, 
imposed by Eq. (28). In itself, the set of Eqs. (27,28) puts no 
constraints on the populations. As the set of Eqs. (27) and (28) 
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always has the trivial solution = 0, any isotropic carrier 
distribution defines at least one stable equilibrium state with 
an associated real polarization defined by Eq. (27). Among 
these equilibrium states, the ground state minimizes the total 
energy and can be determined by inserting Eqs. (27) into the 
variational equation. This yields the condition 



(29) 



relating / and P. Assuming that the interacting ground state 
can be generated adiabatically from the noninteracting tight- 
binding ground state with fk = Pk = 0, the combination of 
the conditional relation and Eq. (27) allows us to express the 
ground state populations in terms of the renormalized energies 
and Rabi energies: 



Pk = 



1 



(30) 



(31) 



Inserting these into Eqs. (22) and (23) yields the coupled set 
of integral equations: 



(ik = \Y.'^{k-k') 



Q.k' 



k' 



(32) 



Sfc = hypk -I'^Vik- k') 



9') with /fc = (1 - sg(Sfc))/2 = 0{kF - k) and E^^ = 
fixes the Fermi wave number kp. Hence, as Efc^o — 0, the 
tight-binding ground state always solves the gap equation and 
hence corresponds to a stationary mean field solution of the 
many body Hamiltonian, though not necessarily the ground 
state. 



B. Properties of the Ground State 

To achieve intuitive insight into the nature of the ground 
state, and to derive a criterion for the existence of a nontrivial 
solution of the gap equation, we implement the conditional 
relation (29) via 



fk — 
Pk^ 



l+/32c/)2(fc) 
l+/3202(fc)2- 



(34) 
(35) 



Without loss of generality, we can assume to be a nor- 
malized wave function and f3 is an additional parameter con- 
trolling the total density. With the aid of the wave function 
one can construct the exciton creation and destruction opera- 
tors 



k' 



B^ = Y,meihik 

k 

B ^^(l){k}h_kek 



(36) 
(37) 



X cos(6i - 61') 1 



(33) 



from which the renormalized and Rabi energies may be cal- 
culated numerically. The set of Eqs. 32 and 33 are equivalent 
to those derived in Ref. 16 and combines the so called gap- 
equation for Clk with the equation for the renormalized single 
particle energy Efc. 

The set of Eqs. (30)-(33) displays several interesting fea- 
tures. From Eq. (31) one recognizes, that fk > 1/2 implies a 
negative renormalized energy < and fk < 1/2 implies 
Sfc > 0, respectively. These relations can also be deduced di- 
rectly from equation 27 by noting that both Pk and Clk are 
positive quantities. From Eqs. (32) and (33), we notice the 
properties Clk=o > ^k and Sfc=o = for any isotropic par- 
ticle distribution, the latter relation resulting from angle inte- 
gration. Hence, for any nontrivial isotropic solution of the gap 
equation one has fk=o = Pk=o ~ 1/2, i.e. the electrons are 
in a state that mixes the tight-binding valence and conduction 
bands with equal probability. 

In general, the solution of the gap equation is not unique. 
This reflects the fact that the many-body Hamiltonian may 
have more than one stationary mean-field solution fulfilling 
the variational condition (26). In particular, the gap equation 
(32) always has the trivial solution ^k = Pk = 0, correspond- 
ing to a completely incoherent state. For the incoherent state, 
Eq. (31) simplifies to = hvpk — '^^f V{k — fe') cos{9 — 



that generate the transformation U{/3) = exp{f3B^). Acting 
on the tight-binding ground state, U creates a coherent exciton 
state 



|*)bcs = Ce^^' |^')tb = l[(uk + vkhl^el) |^) 



/TB, 



(38) 

reproducing the expectation values (34) and (35). Here, 
the normalization constant is C = exp{— 1/2^^ ln(l + 
(j)'^ {k)} . The coherent excitonic state is equivalent to 
the BCS state with Uk = 1/V^ 
/3(l,{k)/./TTW¥{k)- 



f3'^cj)\{k) and Vk 



As any Hartree-Fock state is uniquely determined by the 
expectation value of the basic single-particle operators, Eq. 
(38) represents the ground state if we identify — 
fk I Pk for any nontrivial solution of the gap equation and fix 
/3 by normalizing (\>. At low densities, the commutator 



B.B 



>2(fe) 1 



(39) 



is quasi-Bosonic and the exciton state is formally equivalent 
to a coherent Glauber state in quantum optics. 

Clearly, for any fixed set {(/()(fc)}, the total particle density is 
a monotonously increasing function of /? and the tight-binding 
ground state with Pk = fk = {H)hf — corresponds to 
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FIG. 5. Results from evaluating the energy (a) and the energy per 
particle (b) for the state (38) for different values of /?, with a fixed 
wave function ^(fc). 



/3 = 0. For small 13, the lowest-order contribution to the total 
energy is given by 



E{f3) = 



kk' 



(40) 



Subsequent variation of E/[3'^ — ji 4>{kY yields the linear 
Wannier equation for the wave function with eigenvalue ji. 
In the limit of vanishing density, the system energy is i? = 
— N jjL, where N — (i"^ is, the total number of pairs in 
the low density limit. Hence, if the linear Wannier equation 
allows for bound-state solutions with /i < 0, the system can 
gain energy by forming bound electron-hole-pair states and a 
nontrivial solution of the gap equation is expected to exist. 

In Ref. 30, bound exciton solutions are reported for parame- 
ter conditions coiTesponding to the regime of strong Coulomb 
interaction. These solutions vanish as the system transitions 
into a weakly interacting regime where the graphene fine- 
structure constant ac < 1/2. To demonstrate the properties 
of the excitonic transformation (38), we take the wave func- 
tion 4){k) obtained as numerical solution of the Wannier equa- 
tion, and evaluate the energy of the state (38) for different val- 
ues of (3. For ao = 2.4, which coiTesponds to the coupling 
strength in vacuum for vanishing dynamical screening, the re- 
sult is shown in the left pait of Fig. 5. Increasing /3 increases 
the number of pairs contributing to an energy reduction while 
simultaneously the energy gain per pair (Fig. 5b) decreases 
due to the Pauli blocking. As a result, the obtained energy 
function has a clear minimum at Pcd ~ 0.2. We can also note 
that this minimum is clearly below the energy at /? = 0. The 
energy per pair (Fig. 5b) at f3 — /3c is approximately one half 
of the exciton binding energy (E) / (N) |^^q. 

As any nontrivial solution of the gap equation produces a 
wave function via f3(f>{k) = fk/Pk^ solving the gap equation 
is equivalent to the minimization of the ground state energy 
with respect to the full set f3(j){k) and must yield better results 
than variation with respect to the single variational parameter 
(3. Hence, the existence of bound-state solutions of the linear 
Wannier equation is an unambiguous sign for a nontrivial so- 
lution of the gap equation. 

The solution of the gap equation defines the mean field 
Hamiltonian 



- n. 



k^k 



Skh- 



^1 



0.2 




1 J 








0.1 
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FIG. 6. Numerical solutions of the gap equations, Eqs. (32) and (33) 
for qg = 1A. (a) Ground state Rabi energy Q. (dashed line), ground 
state renormalized energy E (dash-dotted line), and — (Sl^ + 
Y?)^!"^ (solid line), all versus kd. Energies in units of Eq — hvp /d. 
(b) Ground state population / (solid line) and polarization P (dotted 
line) versus kd. These are related to the energies in (a) via Eqs. (30) 
and (31). Note how both / and P go to 1/2 at fc = 0, and that the 
region where / > 1/2 corresponds to the region where E < 0. 



with the single particle dispersion 



± 



(41) 



(42) 



The new valence band is below the tight-binding valence 
band if either ftk — and = has a solution with 
positive Fermi wave number, in which case the new valence 
and conduction bands are degenerated at the Fermi level at 
fc = and k = kp, or the gap equation has a nontrivial so- 
lution with Clk 7^ 0, for which the spectrum exhibits a gap 
of magnitude 2^1^=0 at the Dirac points. Obviously, or rather 
by construction, the BCS state is the ground state of the mean 
field Hamiltonian, and the energy required to add an electron 
in the new conduction or valence band is given by E^'' and 
that to create an electron-hole pair by + EJ,. Neverthe- 
less, the density of states in the pseudogap is nonzero since 
the Bogoliubov state characterized by the order parameters 
{fifejSfc} has a finite overlap to states with order parame- 
ters { fife + Sflk , Sfc + (5Efc } . Hence, within the Bogoliubov 
basis, single-particle excitations are described by scattering 
within the pseudoparticle bands, while collective excitations 
of the many-body system lead to deformation of the bands that 
vary the system energy continuously through the pseudogap. 



C. Numerical Solution of tlie Gap Equations 

Fig. 6 shows the solution of the full gap equations for 
ac — 2.4. Numerically, we evaluated the gap equation itera- 
tively, using the optimized excitonic state generated by the so- 
lution of the linear Wannier equation as initial guess. Conver- 
gence was obtained within less than 10 steps. We checked the 
robustness of the numerical solution against the initial guess. 
Starting with different values of /? requires only a few more 
iterations until convergence is reached but produces the same 
solution. Due to the form factor in the Coulomb matrix ele- 
ments, all integrals are well behaved and no additional cut-off 
parameters are required to achieve convergence. Fig. 6a shows 
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FIG. 7. Solutions of the gap equation for different values of the fine- 
structure constant. Upper part: (a) ground state Rabi energy Q, and (b) 
renormalized single particle energy E versus kd for uq = 0.9 (solid 
yellow line), 1.0 (long-dashed green line), 1.1 (medium-dashed blue 
line), and 1.2 (short-dashed purple line). Energy values have been 
obtained using using vp ~ 9.07 ■ 10'^ and d = 0.18 A. Lower part: 
population (c) and polarization (d) distributions. 



the resulting internal field and the renormalized energies in 
units of Eq — hvp/d. 

We see that the full solution of the gap equation produces a 
large internal field which decreases monotonously as function 
of k. Overall, the renormalized energy deviates only slightly 
from the bare single particle energy, which is due to the angle 
dependence of the Coulomb integrals. At fc = 0, the renormal- 
ized energy starts at zero and for small positive fc, it becomes 
slightly negative. In this region /(fc) > 1/2, as can be seen 
from Fig. 6b showing the population and polarization distri- 
butions. At fc = 0, both 7=1/2 and P = 1/2. Regarding 
the first derivative of the renormalized energy as renormalized 
Fermi velocity, our calculations do not confirm a logarithmic 
divergence predicted in Refs. 2 and 8 but instead, the Fermi 
velocity decreases and even becomes negative. 

In Fig. 7, we show the internal field (Fig. 7a) and renor- 
malized energies (Fig. 7b) for ac — 0.9, 1.0, 1.1, and 1.2. 
From the Wannier equation, we know that the exciton bind- 
ing energy decreases rapidly with decreasing effective fine- 
structure constant, and from this, we expect a similar behavior 
for the internal field. Indeed, decreasing the coupling constant 
from 1.2 to 0.9 decreases the internal field about one order 
of magnitude, while the renormalized energies are hardly dis- 
tinguishable from the bare single-particle energies. Under all 
conditions, the particle and polarization distribution at fc = 
are 1/2 see Figs. 7 c and d. With increasing wave number, the 
distributions fall off very rapidly and the integrated particle 
and polarization density decrease with decreasing ac- 

In Fig. 8, we show the total energy density (Fig. 8a), en- 
ergy gain per pair (Fig. 8b) and pair density (Fig. 8c) as 
function of ac- The energy gain per pair starts at zero for 
ac — 1/2 and increases rapidly with increasing ac, yield- 
ing approximately 1.6 eV for ac = 2.4. Increasing ac from 
.5 to 2.4 increases the total energy gain per cm^ over ap- 
proximately 20 orders of magnitude, which is a combined ef- 




FIG. 8. Energy density (a), energy gain per pair (b), and pair density 
(c) as function of the coupling constant aa- The dotted curves show 
the quantities on a logarithmic scale (left axis) and the solid curves 
present the same quantities on a linear scale (right axis). The stars 
correspond to the gap equations solutions in Fig. 7. 



feet of the increasing energy gain per pair and the increas- 
ing pair density. At ac = 2.41 the pair density is (N) j A = 
2.16 • lO^'^cm^^, which should be compared with the total 
density uq = 1.27 




-2 





ky- Ky{l/nm) 



FIG. 9. Dependence of the quasiparticle band structure El. — El. 
(42) on the effective fine-structure constant aa- The curves shown 
are for ac = 0.9 (solid yellow line), 1.0 (long-dashed, green line), 
1.1 (medium-dashed, blue line), and 1.2 (short-dashed, purple line), 
together with the linear tight-binding dispersion (dotted black line). 
The shaded, dotted area shows the full tight-binding dispersion. 
These band structures correspond to the gap equations solutions in 
Fig. 7. 



Fig. 9 shows the resulting quasiparticle dispersion for dif- 
ferent values of the effective fine-structure constant in com- 
parison with the noninteracting TB bands. For values of ac 
exceeding the critical value ac = 1/2, a gap opens that in- 
creases with with increasing ac- For values of ac < 1.1, the 
Bogoliubov bands lie well within the validity range of the lin- 
ear approximation to the tight-binding bands. For higher val- 
ues, corrections beyond the linear approximation should be 
taken into account, however, this is beyond the scope of this 
article. 
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IV. OPTICAL RESPONSE 

To illustrate the influence of Coulomb effects on mea- 
surable quantities, we calculate the optical response of our 
graphene model system. For this purpose, we derive the 
graphene Bloch equations (GBE) and compute the linear re- 
sponse to an external optical field. In all calculations, the ini- 
tial conditions will be chosen such that the graphene sheet is 
excited from the BCS ground state by a normally incident op- 
tical or terahertz pulse. 

As shown in section 11, at normal incidence, the two de- 
generate Dirac points couple to the right and left circular po- 
larized components of the transverse optical field. Hence, the 
circular polarization components provide a natural basis for 
the optical field in which the different polarizations are de- 
coupled. Each polarization component of the vector potential 
obeys the wave equation 



C2 dt^ 



in 



(43) 



Here, j- 
current 



c{ g^±t ) is the expectation value of the total 



niQC 
e\TT 
mo 



-^A±(z^O,t)/(z)+j±/(z) 



(44) 



with f{z) — J (Pp\(t>p^ [f)]"^ ~ on the scale of the opti- 
cal or terahertz wavelength. The total current consists of two 
contributions. The first is simply proportional to the vector po- 
tential and hence purely classical. Its only dependence on ma- 
terial quantities is via the charge density p{r) = —enof{z), 
where no is the total density of the electrons in the 7r-band. 
This term is responsible for a Drude-like response. The second 
part of the current is purely particle like and associated with 
inter- and intraband transitions. Due to the phase 6, isotropic 
electron and hole distributions or polarizations do not con- 
tribute to a macroscopic particle current. Hence, if the system 
is in its ground state and in the absence of external fields, there 
is no macroscopic current. 



A. Graphene Bloch Equations 

To calculate the particle current, we divide the polarization 
and the particle distributions into a static part arising from the 
ground state, and a dynamical part arising from the optical 
excitations; 



fk = fk + ^fk 



(45) 
(46) 



and correspondingly for the renormalized and Rabi energies. 
In a first approximation, all contributions beyond Hartree- 
Fock are treated phenomenologically by introducing constant 



dephasing and relaxation rates, i.e. ^AP^ = — 72APfc and 
similar for the populations. Noting that P^ is real, the GBE 
for the dynamical part of the polarization and particle distri- 
bution are obtained as 



ih-APk = 2 (Sfc + AE^) APfc - (1 - 2/fc - 2A/fc)Ar!f 

+ 2ASfc (Pfc + APfc) + z(l - 2/fc - 2A/fc)Al]fc 
+ 2PfcAE^ + 2f^fcA/fc - i%2APfc, (47) 

+ 2Ar!fc (Pfc + Re [APfc]) 

- 21m [PkAQ^ + QkAP^] - h^nh. (48) 

In Eq. (47), the first line is the direct generalization of the 
familiar homogeneous part of the SBE, producing the gener- 
alized Wannier equation. The dynamical part of the Rabi en- 
ergy Ail§ contains the optical field and the internal field of 
the optically induced interband polarization only, while both 
the ground state and the dynamical populations contribute to 
the renormalization of the single particle energies and the 
phase space filling, reducing the effective Coulomb interac- 
tion. The second line describes an energy renormalization due 
to the optical field and Auger scattering proportional to AE^, 
and Auger contributions proportional to Afife that act as a 
source/drain for the dynamical polarization. These processes 
are very ineffective in a wide gap semiconductor and are usu- 
ally neglected in the SBE. The last line Eq. (47) describes the 
coupling of the polarization to the populations via the ground 
state polarizations. These contributions do not exist in the SBE 
and are specific for an excitonic ground state. 

The first line in the equation of motion for the populations, 
Eq. (48), is the direct generalization of the SBE. The respec- 
tive terms contain dynamical quantities only and are at least of 
second order in the optical field. The second line describes the 
conversion of polarizations into populations via Auger scat- 
tering processes proportional to Afl^. The last line describes 
ground state polarization assisted sources, specific for the ex- 
citonic ground state. 

The major effect of the ground state polarization is the in- 
troduction of a linear source for the optically induced pop- 
ulations /fc. As a result, within a power expansion in terms 
of the exciting field, both the dynamical polarization and the 
population contain all powers of the exciting field. For any 
arbitrary order, the equations of motion are inherently cou- 
pled by the contributions proportional to P^ in the last lines 
of Eqs. (47) and (48), and hence, must be solved simultane- 
ously rather than iteratively. It is exactly this polarization me- 
diated coupling that is responsible for the occurrence of the 
Bogoliubov gap in the linear spectrum. 



B. Linear Optical Spectra 

To illustrate the basic effects of several contributions to the 
optical response, we solve the linearized GBE for a given ex- 
ternal optical field. Examples of the results for ac — 2.4 and 
h-f — 5 ■ 10^'*£'o are shown in Fig. 10. The spectrum shows 
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huj/Eo 

FIG. 10. Imaginary part of the linear susceptibility for a graphene- 
like system with ao = 2.4. Energy in units of Eo = fivp/d. 

ATTlm[j^ /ojA^] which is a direct measure for the absorp- 
tion. We notice pronounced excitonic resonances at low ener- 
gies followed by a spectrally flat response. Due to the optical 
selection rules, all bright excitonic resonances have a p-like 
symmetry. Similar to the spectra in semiconductors, the peak 
height increases with increasing binding energy, showing that 
oscillator strength is transferred to strongly bound excitons. 

Ignoring all Coulomb effects in the GBE, other than those 
responsible for the excitonic ground state, gives the artifi- 
cial spectrum shown by the dashed (blue) curve in Fig. 10. 
This solution is equivalent to that obtained using only the 
mean-field Hamiltonian (41). The resulting spectrum is pro- 
portional to the density of states (DOS) (divided by uj) of the 
Bogoliubov bands and has no additional structure. The onset 
of the continuum absorption starts with a Xj^/huj — i?g-like 
damped singularity exactly at the Bogoliubov gap. The origin 
of the singularity is the shift of the band minimum to finite fc- 
values (see Figs. 6 and 7). Within a quadratic approximation 
E(k) =Eg+ /i^(|fc| - fco)V2mcff, the DOS 

consists of a step-like contribution typical for 2D parabolic 
bands and a singular part proportional to the shift fco that is 
similar to the free-particle result of a semiconductor quantum 
wire.^** Energetically above the Bogoliubov gap, the optical 
response assumes the 'universal' value given by ^TTa{e* / e)"^ 
where a = /he is the fine-structure constant and e* = 
e\-K\ / \/2niQV F is the effective charge, producing the effec- 
tive medium Hamiltonian Hq + Hj = vpcr ■ {p — e* /cA) 
where cr = {ax, ^cTy) combines the Pauli spin matrices. The 



1 

i 


1 

I 




1 





0.1 0.2 

FIG. 1 1 . Effect of dephasing on the spectra of a graphene-like system 
with QfG = 2.4. The curves show the results for dephasing constants 
(from most to least peaked) h'y/Eo = 5 • 10"'', 1 • 10"^, 2 • 10"^, 
5 • 10"^, and 5 • 10~^. Energies in units of Eo = hvp/d. 



effective charge depends on the optical matrix element and 
can either be calculated from the carbonic wave functions or 
be used as a fit parameter. An additional factor 2 arises if the 
spin degeneracy is taken into account. Note that there is no 
valley degeneracy for the excitation with circularly polarized 
light. 

Including all Coulomb interaction terms, we obtain several 
clearly recognizable excitonic resonances below the pseudo- 
gap, shown by the solid (pink) curve in Fig. 10. The reso- 
nances arise from the Coulomb interaction of the optically in- 
duced excitations only, while the excitonic ground state pop- 
ulations open the required gap. 

To study the influence of dephasing, we repeated the cal- 
culation for the dephasing values in the range from hj — 
5 • 10"^i?o (lightest, most peaked curve) to ~ 5 ■ 10~^i?o 
(darkest, flattest curve). As we can see in Fig. 11 the exci- 
tonic spectra broaden and the fine structure close to the gap 
smoothens to a flat, continuum-like response. The lowest ex- 
citon resonance is clearly recognizable for dephasing rates 
up to 5 • IG"^^'^ 

To study the influence of the Auger contributions, we 
present in Fig. 12 results where we artificially vary the rel- 
ative strength of the Auger terms between (no Auger terms) 
and 1 (full Auger contributions). The comparison shows that 
the Auger contributions not only increase the exciton binding 
energy significantly, but they are also responsible for the split- 
ting of the individual resonances. In the absence of the Auger 
contributions, the two distinct projection states of the angu- 
lar momentum onto the propagation direction of the incident 
light are degenerated. As the Auger contributions invert the 
rotational symmetry they remove this degeneracy. The split- 
ting of the resonances increases with increasing strength of 
the Auger terms and is symmetric with respect to the energy 
shift. The oscillator strength is distributed nonequally among 
the resonances. The curves in Fig. 12 show Im [j/^], which 
differs by a factor 1 /u from the absorption spectra. This pre- 
sentation shows that the the heights of the absorption peaks 
vary exactly like l/w, with a constant prefactor that is consid- 
erably larger for the lowest lying state and independent of the 
strength of the Auger terms. As the Auger terms vanish for 
s-like states, the splitting confirms the p-like symmetry of the 
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FIG. 12. Imaginary part of the linear susceptibility for a graphene- 
like system with aa = 2.4. Energy in units of Eq = hvp/d. 




FIG. 13. Imaginary part of the linear susceptibility for a graphene- 
like system for aa = 0.9 (solid yellow line), aa ~ 1.0 (long-dashed 
green line), qq ~ 1.1 (medium-dashed blue line) and qq ~ 1.2 
(short-dashed purple line). These spectra correspond to the gap equa- 
tions solutions in Fig. 7 and to the band structures in Fig. 9. 



optically bright states. 

In Fig. 13 the absorption spectra for several values of the 
coupling strength and a fixed dephasing rate hj = 13 meV 
are shown. Variation of the effective coupling strength can be 
realized by embedding the graphene sheet in, or putting it on 
top of a dielectric medium, altering the static screening e = 
or e = (l+Er) /2 respectively. In the strong Coulombic regime 
with ac > 1/2, all spectra show one clearly recognizable dis- 
crete absorption peak and a flat continuum response. Both the 
spectral positions and the height of the absorption peak de- 
pend extremely sensitive on the effective coupling strength, 
i.e. on screening by the dielectric environment. Qualitatively, 
any single spectrum is similar to that of a conventional semi- 
conductor. However, the dependence on the coupling strength 
is very different. If one compares e.g. quantum well systems 
with a different background screening and defines the effec- 



tive coupling strength as = e^/ehc = a/e, the spectral 
position of the gap is independent of the background refractive 
index while the exciton binding energy varies like a^gm,-c^ 
and mi is the reduced mass of the electron-hole pair. The net 
result is a red shift of the excitonic resonances with increasing 
coupling strength. In graphene, the exciton binding energy, 
defined with respect to the Bogoliubov band edge, increases 
also with increasing coupling strength, but simultaneously, the 
band edge shifts towards higher energies, resulting in a net 
blue shift of the lowest excitonic resonance with increasing 
coupling strength. 

V. SUMMARY AND CONCLUSIONS 

In conclusion, we presented a framework to determine the 
ground state and optical response of graphene and graphene- 
like systems. Our method is based on the equations of motion 
for the basic variables, combined with a variational Ansatz for 
the ground state. Even though we have only presented results 
in linear approximation, our method allows for a selfconsis- 
tent inclusion of e.g. higher order many-body correlations or 
dynamical screening, treating the ground state properties and 
excitation dynamics on equal footing. It can easily be general- 
ized to finite temperature, be adapted to other strongly corre- 
lated systems or to include additional constraints like doping. 

Within the Hartree-Fock approximation, our procedure 
produces the gap equation for the ground state. We could 
relate the criterion for exciton condensation of the ground 
state with the existence of bound exciton solutions of the lin- 
ear Wannier equation discussed in Ref. 30 and solved the 
gap equations numerically. Our results for the ground state 
are in general agreement with both Monte Carlo simulations 
and variational approaches studying the semimetal to insula- 
tor transition '^'"^ and predict a a critical value of a = 0.5 for 
the onset of exciton condensation. 

Consistently with the level of approximation for the ground 
state, we derived the graphene Bloch equations for the 
strongly interacting regime on the singlet level, which allows 
us to calculate the optical response from the BCS ground state. 
Unlike in direct-gap semiconductors, where electron-hole re- 
combinations are dipole allowed, s-like excitons in graphene 
are optically dark and hence radiatively stable. Signatures of 
the excitonic ground state are the existence of a pseudo gap 
and several bright, p-like excitonic resonances in the optical 
spectra, with a spectral position depending extremely sensi- 
tively on the effects of screening. 

Topics of future research are e.g. the inclusion of dynam- 
ical screening, that might reduce the nominal effective fine- 
structure constant ac = 2.4'"* significantly, the inclusion of 
Coulomb scattering, correlations beyond the singlet level, as 
well as optical nonlinearities. 
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